LIBRARY  OF  THE 

UNIVERSITY  OF  ILLINOIS 

AT  URBANA-CHAMPAICN 

510.84 

I&63c 

HO.&1-70- 


AUG.    51976 

Ihe  person  charging  this  material  is  re- 
sponsible for  its  return  to  the  library  from 
which  it  was  withdrawn  on  or  before  the 
Latest  Date  stamped  below. 

Theft,  mutilation,  and  underlining  of  books 
are  reasons  for  disciplinary  action  and  may 
result  in  dismissal  from  the  University. 

UNIVERSITY    OF     ILLINOIS     LIBRARY    AT    URBANA-CHAMPAIGN 


EN 


tm 


jaw.  '  i 
APR     8  I 


fi .'•■  .<! 


n,r 


torn 


PHOTO  REPRODICTION 

OCT  221ECD 

PHOTO  REPRODUCTION 

NOV  1  6  ftflf 

reproduce* 


ft 


OCT 

EHOIO 


8S0I0 

OCT 


: 


o  i  ffn 


f;?fl 


2  7  KC1 

^PRODUCTION 

TjRODUtf: 

3  RECTI 


L161  —  O-1096 


Digitized  by  the  Internet  Archive 

in  2012  with  funding  from 

University  of  Illinois  Urbana-Champaign 


http://archive.org/details/simultaneousfitt61belf 


ENGINEERING  LiBRARV 
UNIVERSITY  OF  ILLINOIS 


meed  Computation 


UNIVERSITY  OF  ILLINOIS  AT  URBANA-CHAMPAIGN 

URBANA,  ILLINOIS  61801 


CAC  Document  No.  6l 

SIMULTANEOUS  FITTING  OF  EXPONENTIAL 
DECAY  CURVES 


By 


Geneva  G.  Belford 


April  1,  1973 


-2' 


CAC  Document  No.  6l 


SIMULTANEOUS  FITTING  OF  EXPONENTIAL 
DECAY  CURVES 


by 

Geneva  G.    Bel ford 


Center  for  Advanced  Computation 

University  of  Illinois   at   Urbana-Champaign 

Urbana,    Illinois      6l801 


April  1,    1973 


This  work  was    supported  in  part  by  the  Advanced  Research  Projects 
Agency  of  the  Department  of  Defense  and  was  monitored  by  the  U.    S, 
Army  Research  Office-Durham  under  Contract   No.    DAHC0U-72-C-0001. 


bio, 84-  0 

ABSTRACT:    This  paper  deals  with  characterization  of  best  approxima- 
tions to  vector-valued  functions.   The  approximations  are  themselves 
vector- valued  functions  with  components  depending  nonlinearly  on  the 
approximation  parameters.   The  constraint  is  imposed  that  certain  of 
the  parameters  should  be  identical  for  all  components.   An  applica- 
tion to  exponential  approximation  is  discussed  in  some  detail. 
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1.   Introduction 

The  work  reported  in  this  paper  was  motivated  by  the  following 
problem:   Suppose  a  set  of  experimentally  determined  exponential  decay 
curves  is  given.   It  is  desired  to  approximate  the  curves  by  functions 
of  the  form  a  exp  (8x),  where  3  should  be  the  same  for  the  entire  set  of 
curves  and  a  may  vary  from  curve  to  curve.   The  problem  is  to  determine 
how  such  approximation  might  best  be  made.   This  problem  arises  in  a 
number  of  physical  situations.   In  chemical  kinetics,  for  example, 
monitoring  of  a  chemical  reaction  which  obeys  a  first-order  rate  law 
leads  to  just  such  exponential  data,  from  which  one  wishes  to  extract 
a  best  8  although  the  initial  amount  of  material  (a)  varies  from  experi- 
ment to  experiment . 

In  previous  papers  [l,  2,  3]  this  type  of  constrained  vector- 
valued  approximation  was  studied  for  the  simpler  situation  where  the 
approximating  functions  depend  linearly  on  the  parameters.   In  this  paper, 
results  for  nonlinear  approximation  are  presented.   Section  2  contains 
a  precise  formulation  of  the  problem  and  characterization  theorems 
applicable  to  the  construction  of  best  approximations  from  general 
classes  of  nonlinear  families.   In  §3  the  particular  problem  posed 
in  the  preceding  paragraph  is  discussed  and  a  very  simple  alternation 
theorem  is  obtained.   Finally,  §k    describes  how  the  theoretical  results 
have  been  used  to  design  a  successful  computer  algorithm  for  the  simul- 
taneous fitting  of  exponential  decay  curves. 


2.      Formulation  of  the   Problem  and  Characterization  Theorems 

Let   g    ,    g   , ...    ,    g     be  a  given  set   of  real   functions   continuous 

on  a  closed  interval   I  of  the  real  line   and  let   g  denote  the  ^-dimensional 

vector-valued  function  with   components   {g.}.      Let  V  be  an  n-paarameter 

family  of  functions   in  C(l).      Denote  an  arbitrary   element   of  V  by 

<j>(a,,    •••>   a    ;    x),   where    (a,,    ...,   a    )    e   R     is  the   parameter  vector, 
Y     1  n  1  n 

and  assume  that   $   depends   continuously  on  the  a.    as  well   as  on  x.      For 
any   integer  m   (0  <_  m  <_  n),   we  then   define  the   family  of  approximating 
vectors 

F  =   {(f1(a,x),f2(a,x),...,    f      (a,x)    ): 

a=    (air   a12,    ....    alm,a2r    ....    a^,    ....   a^.a^,    ...,-c.J 

eR   ,    q  =  m(£-l)+n,    and 

f .  (a,x)   =   <b(a.,,    •••»   a.    ,   a    , ,  ,    ...,   a    ;x)} 
l  Y     ll  lm       m+1  n 

For  example,   the  approximating  family  for  a  set   of  exponential   decay 

curves    (as   described  in  the   introduction)   would  be   defined  by 

4>  =   a-,exp(apx)    and  m  =  1.      In  order  to   avoid  the   double-subscript 

notation,   we  will  henceforth  write  a  =    (a, ,ap,    ...,    a    ). 

The  norm  used  in  this   paper   is  the  usual   uniform  norm;    that 

is,    if  f  is   any  vector- valued  function  with   components    f.    in   C(l),   the 

norm  N(f)    is   defined  by 

N(f)   =  max     ||f.  || 
i  1 

where  ||f.||=max    |f.(x)|    . 

1  xeI        ■L 

An  element   f  in  F  is  then   called  a  best   approximation  to   g  from  F  if 

N(g-f)   =  inf  N(g-f)    E  p(g)    . 
feP 


One   does  not    in  general   expect   the  existence  of  such   a  best 
approximation   since  best   approximations   to  a  single   function   from  non- 
linear  families  often   fail  to  exist.      Uniqueness  of  the  best   approxima- 
tion  is   also  the  exception   rather  than  the  rule,    as   is   the   case   for  the 
simpler  situation  when  F   is   a  linear  subspace    [l].      Comments   on   existence 
and  uniqueness   can  be  provided  in  particular  cases,   however,    as  will  be 
seen   in  the  next    section. 

Our   first   theorem,    analogous  to   Theorem  85  of  Meinardus    [8]    and 
having  very  similar  proof,    provides   a  lower  bound  on   p(g). 

Theorem  1:      Suppose  that    for  some  a  e   R     and  subsets 
DkC  I    (k  =  1,    2,    ...,    l)y    f   (a,x)   -  gfcU)   4  0   for  all   x  eD     and  all   k. 

Suppose  also  that  there   exists   no  parameter  b   e  R      such  that 

(1)  [gk(x)   -   fk(a,x)]    [fk(b,x)   -    fk(a,x)]   >    0   for  all   x  e   Dfc 
and  all  k. 

Ihen 

(2)  p(g)   >   inf  inf      |g,(x)   -   f,(a,x)|. 

k     xeDk       K  K 

Specializing  the  sets  D    ,    one   is   then  led   immediately  to   a 
sufficient   condition   for  a  best   approximation. 

Definition  1:      The   pair    (x,k)    is    called  an   extremal   of  the 
approximation   f (a, . )   to  g   if 

|gk(x)   -    fk(a,x)|      =  N(g  -   f(a,    .)    ). 
Now  let   Ej(a)      =    {x   :    (x,k)    is   an  extremal}    . 

Theorem  2:      If  there   exists  no  b  e   Rq  such  that    (l)   holds    for 
all   x  s   ^v(a)    an^  all  k,    then   f(a,x)    is   a  best   approximation  to   g  from  F. 


h 


In  order  to  obtain  a  necessary  condition  for  a  best  approxi- 
mation, additional  hypotheses  on  F  are  needed.   Henceforth,  we  assume 

that  the  partial  derivatives   9f,  ,    s  all  exist  and  are  continuous  for 

k(a,x) 

9a. 
J 

(a,x)  e  R   x  I.   The  following  theorem  then  holds.   Again,  this  is  a 
straightforward  vector  analog  of  a  known  theorem  [8,  p.  lUO-l],  and 
the  proof  will  be  omitted. 

Theorem  3:   Let  f(a,x)  be  a  best  approximation  to  g  and  assume 
that  p(g)  >  0.   Then  there  exists  no  b  =  ($  ,  ...,  3  )e  R  such  that  for 


all  k 


(3) 


gk(x)  -  fk(a,x) 


I   6,  3fk(a>x) 


j=l  "    9a 
J 


>  0  for  all  x  e  E  (a). 


Notice  that  both  the  sufficient   and  the  necessary   conditions 
for  a  best   approximation   involve   systems   of  inequalities    ((l)    and   (3), 
respectively).      Computationally,   the  linear  inequalities    (3)   are  much 
easier  to  handle  than  the  nonlinear  ones    (l).      It   is  then  of  interest 
to   determine  whether  under  some   conditions   the  linear  inequalities    (3) 
may  also   figure   in  a   sufficiency  theorem.      For  approximation  of  a  single 
function,    Krabs    [6,7]   has   discussed  several   such   conditions,   one  of  which, 
the  representation  condition    [6],    readily   extends   to  the   vectorial   case 
and  yields   results   applicable  to  the   exponential   approximation  problem 
in  which  we   are   interested. 

Definition  2:      The    family  F  is    said  to   satisfy  the  representation 

# 

condition   if  for  every  pair  of  functions    f(a,x),    f(b,x)    in  F  there  exist 

real  numbers    c.(a,b)    (j=l,    ...,    q)    and  functions    <J>    (a,b;x)    (k  =  1,    ...,    l) 

3  k 

positive  on   I  such  that 


q        .   , 

(k)      f   (a,x)  -  f(b,x)  =  ♦.(a,b;ac)  I  c.(a,b)  8fkU'x)   for 
k  K         K       ,1=1  J        3a. 

k  =  1,  . .. ,  £.  J 

The  required  connection  between  inequalities  (l )  and  (3) 
is  then  provided  by  the  following  theorem,  which  is  an  extension  of 
Krabs'  Theorem  3  [6]  and  has  a  virtually  identical  proof. 

Theorem  h:      Let  F  satisfy  the  representation  condition  and 
suppose  that  for  some  set  of  closed  subsets  D.  £  I  (k  =  1,  2,  .  ..,  i) 

there  exists  no  b  =  (3_,  ....  $  )  such  that  (3)  holds  for  all  x  in  D, 

1  q  k. 

and  all   k.      Then   for  all  b   e  Rq, 

(5)      min  min    [g.  (x)   -    f.  (a,x)]    [f,  (b,x)   -    f   (a,x)]    <_  0. 
k     xeLV     k  k  k  K 

IV 

The  conclusion  of  this  theorem  is,  of  course,  equivalent  to 
the  statement  that  for  no  b  does  (l)  hold  for  all  x  in  D  and  all  k. 

Putting  together  the  preceding  results,  we  obtain  the  follow- 
ing useful  characterization  theorem. 

Theorem  5-      Let  F  satisfy  the  representation  condition  and 
suppose  that  p(g)  >  0.   Then  f(a,x)  is  a  best  approximation  to  g  if 
ana  only  if  there  exists  no  b  =  (^  ...,  Bq)  e^  such  that  for  all  k 
and  all  x  e  E  (a)  the  inequality  (3)  holds. 

In  applying  this  characterization  theorem,  it  is  useful  to 
keep  two  things  in  mind.   Firstly,  consistency  or  inconsistency  of  the 
inequalities  (3)  does  not  depend  on  the  magnitude  of  the  approximation 
error  (N  (g  -  f(a,  .)))  but  only  on  the  signs  sgn[g  (x)  -  f  (a,x)],  which 

K.  K 

we  shall  denote  by  o ,  (x).   Secondly,  a  necessary  and  sufficient  condition 

K. 

for  the   inconsistency  of  the   set   of  inequalities    (3)    (where    (x,k)   runs 
over  all   extremals)    is   that   the  origin    6  =(0,    ...,    0)   of  R     should  lie   in 
the   convex  hull   of  the   set   of  q-vectors 


9fk(a,x)                   3fk(a,x) 
(6)      S  =   {      a    (x)    ( 5 ,    ...,   ~ )    :    (x,k)    is   an   extremal) 

K  oOt  dot 

1  q. 

(This   result   on  linear   inequalities  may  be   found  in  Cheney's  book   [5, 
p.    19].)      This    "convex  hull"   condition  yields   a  set   of  linear,   homoge- 
neous  algebraic   equations,    from  which  one  may  often  quickly  determine 
whether  or  not   a  given  approximation   is  best.      These   equations  may  also 
be  used  to  obtain  alternation  theorems    [5,   p.    75]    a  fact  which  will  be 
exploited  in  the  next   section. 


3.      Application  to  Constrained  Exponential  Approximation 

In  this   section  we   consider  the  particular  case  of  approxi- 

2 
mation  "by  one-term  exponentials    (elements  of  E     =   {aexp(Bx)    :    (a,B)£R   }), 

"with  the   exponential    factor   B   constrained  to  be  the   same   for  all   com- 
ponents.     That    is,    as   noted  before,   we  take   f     =  a,exp(a   x)    for  k  =  1, 
...,£  =   q-1;    the   resulting  family  of  vector-valued  approximating  functions 

will  be   denoted  by   F        .      For  simplicity,   we  shall  take  the   interval   I 

exp 

to  be  [0,  l].   The  existence  of  best  approximations  from  F    is  readily 

demonstrated.   The  compactness  of  any  bounded  set  {  f  :  feF    ;  N(f ) <  M  } 

exp      — 

is  easily  deduced  from  the  known  compactness  result  for  I   =   1  [9].   The 
usual  existence  argument  then  goes  through. 

In  order  to  apply  all  of  the  theorems  of  the  last  section,  in 
particular  Theorem  5,  we  must  first  verify  that  the  representation  con- 
dition holds  for  F    .   (The  smoothness  condition  prefacing  Theorem  3 

exp 

is   clearly  satisfied.)      The  representation   condition   is   known  to  hold 
for  exponential   functions[6],   but   extension   from  scalar  to  vector-valued 
functions   is   certainly  not  obvious,   because   of  the   requirement  that  the 
coefficients    c     not   vary  with  k.      Letting  a  =    (a    ,    ...,    a    )    and  b  =    (6,, 
...,    B    ),   we  need   (from  (h)) 

(7)   a  eXp(a  x)    -   B  exp(B  x)   =   exp(a  x)    4>   (a,b;x)    {c    (a,b)   + 

K.  Q  JC  (^  (^  K  K. 

a  c    (a,b)x   }   for  k  =  1,2,    ...»    q-1. 

K   q 

For   4>     to  be  nonvanishing  on   [0,1  ]   as   required,   the  linear  factor  on  the 
right    (in  braces)  must  have  a  zero  at  the   same  point   as   does  the  left 
side.      This   condition  leads  to 

(8)  \  i  \  -  -\lo«Vek)/<6q  -  v 

(k  =  1,2 q  -!)• 


(Equation    (8)    is   obtained  under  the  assumptions   that   a   6     >   0    (k  =  1, 
•  ••,    q  -  l)    and   3     -   o.     ^  0.      The  other   cases    are  readily  handled  by 
similar  arguments.  )      Furthermore,   the  posit ivity  of  <j>     requires   also 

K. 

that   the   signs  of  both   sides   of   (7)   should  match  at   any  point  x   . 
Taking  x  =  0,   ve  arrive  at  the   condition 

(9)      sgn    (c^  -   3k)   =  sgn    (cfe)      (k  =  1,    ...,    q  -1 ) . 

Now  since   sgn    (a     -   3,  )    =   sgn   [a  log(a ,/3,  )],    it    is   clear  that  by  choos- 

K.  K.  JtC  X       K. 

ing  any  c   such  that  sgn  c  =-sgn(3  -a)  and  then  solving  (8)  for 

Si  M.  St  St 

c    ,    .  ..,    c        »    a  set   of  scalars   c    ,    .  ..,    c      satisfying   (8)   and    (9)   may 
be   found.      The   functions    <f>      (a,b)    are  then   defined  to  be 

K. 

a     exp(a  x)   -    3,     exp(3  x) 
\     ~       exp(aqx)    {ck  +  Vqx}        ,   k  =  1 ,    . .  . ,    q  -  1, 

and  the   representation   condition   is   verified. 

Then   in  order  to  obtain   an   alternation  theorem  characterizing 
best   approximations,   we  may  apply  the   "convex  hull"   condition  as   discussed 
at   the   end  of    §2.      First   note  that  by  Caratheodory 's   Theorem   [5,   p.    IT], 
the   condition  that  the  origin      6   of  R     should  be   in  the  convex  hull   of 
S  may  be  replaced  by  the   condition  that     6    should  be   a  convex  linear 
combination  of  some   q  +  1    (or  fewer)   elements   of  S.      For  our  exponential 
approximation,   then,    the   condition   is   that  there   exist   extremals 
(x    .,k)    (with  k  =  l,    ...,    I,    i   =  l,    ...,    v    ,    and  J  v,     <      q  +  l) 
and  nonnegative   constants   X        satisfying 

K.  1 

J  J   X        =1   such  that 
..     .      ki 
k  l 

Vk 

(10)         ^      Aki      °ki   eXp(aqXki)    =   °    (k  =  1,    ..-,    q   -  1) 
i=l  4 


1-1  \ 

I        1      Aki  aki  ak  \±   eXP(aqXki)  =  °  * 
k=l  i=l 

(Here  we  have  used  a  .  to  denote  a  (x  . ). ) 

ki  K-  ki 

Of  course,  only  those  extremals  (x,->  k)  for  which  A   is  non-zero  play 
a  role  in  actually  characterizing  a  best  approximation.   Thus  one  imme- 
diately sees  from  (10)  that  any  k  for  which  v  =  1  does  not  enter  into 
the  characterization.   Considering  the  various  possibilities  involving 
indices   k  for  which  v  >  1,  one  quickly  arrives  at  the  following 

K. 

alternation  theorem. 

Theorem  6:   The  vector-valued  function  f  is  a  best  approxima- 
tion from  F    to  g  on  [0,1  ]  if  and  only  if  one  of  the  following  con- 
ditions holds. 

(i)  For  some  index  k,  f  is  a  best  unconstrained  approximation 

K. 

to   g     from  E   •    i.e.,   there  are  three  points   of  alternation   if  a,    ^  0   and 
two  points   of  alternation   if  a     =  0   [8,   p.    178]    and     ||  gk  -  f-k\\=  N(g-f). 

(II)      There   exist   two   indices    (say  k  =  1,2)  with   four  associated 
extremals    (x^l),    (x12,l),    (x21,2),    (x22,2)   such  that   0^  ?  0,    c*2  ?  0 

Xll  <    X12'    X21  <   X22'    and 

°11   °12  =  -1 

°21   °22   =   -1 

°11    CT21   =   ~Sgn    ^ai   a2^' 

Example:      Let   g  =    (l,x).      The  best   approximation   from  F  to 

g  on    [0,1]    is   given  by   f     =  -  e        ,    f     =  -  e   gX  with   6  =  log  2.      One 
readily  verifies  that   N(f-g)   =  1/3,    extremals   are    (0,1 ),    (l,l),    (0,2), 
(l,2),    and  the  alternation  requirements   of  condition    (il)   of  Theorem  6 
are   satisfied. 


10 


Finally,  we  take  up  the  question  of  uniqueness.   It  is  obvious 

that  the  factors  a   (k  =  1,  2,  . ..,  q  -  l)  are  not  in  general  unique. 

However,  it  is  likely  to  he  the  parameter  a  that  is  of  principal  interest, 

and,  as  the  following  theorem  shows,  this  parameter  is  (with  a  trivial 

exception)  uniquely  determined. 

Theorem  J:   Let  f(a,x)  be  a  best  approximation  from  F    to  g. 

exp 

Let  a  =  (a  ,  ...,  a  ).   Suppose  that  either  (i)  condition  (i)  of  Theorem 

6  holds  for  some  k  such  that  a  4   0,  or  (ii)  condition  (il)  of  Theorem  6 

holds  for  some  pair  of  indices  k  ,  k  such  that  a   ^  0,  a   ^  0.   Then 

1   d  12 

if  f(b,x),  with  b  =  (3  ,  ...»  B  ),  is  any  other  best  approximation, 
a  =  6  . 

q.   q. 

The  proof  of  this  theorem  may  be   found  in  detail   in    [h] . 


11 


h.      Construction  of  Best   Simultaneous   Exponential  Approximations 

The  problem  of  finding  a  best    approximation  to   a  general   g 

is    simplified  enormously  by  the  knowledge    (from  Theorem  6)   that   one 

need  not   consider  more  than  two   of  the   component   functions   simultaneously. 

Notice  that  the  key  problem  is  to   determine  a  best  a    ,    since  then  a  set 

of  a,  's    (k<q)   may  be   determined  from  the   condition  that  a,     exp(a   x) 
k  k  q 

should  be  the  best   approximation  to   g     from  the  linear   family 

{aexp(a   x)    :    aeR}.      Thus   one  would  proceed  by   constructing  best 

unconstrained  approximations   from  E     to   each  g .      If  condition   (i)   of 

Theorem  6  is  to  hold,    the  best   exponential   factor  must  be  that   for  the 

best   unconstrained  approximation  f     such  that      |f     -  g    | |    =  max      |f     - 

g    I  I .      One  then  checks   to  see  whether,    using  this   exponential   factor  as 
J 

a        the  other  curves  may  be   fitted  to  within   error      If,     -   g,  I  I  .      If 
q  J  '  '    k        TcM 

not,    the  next   step   is  to   examine   all   pairs   g,     ,    g,     •      For  many  of  these 

pairs,   the  best   approximation  is   characterized  by   condition    (I)   of 

Theorem  6  and  is  therefore  of  no   further  interest.      If  best   approximations 

f     ,    f        characterized  by   condition    (il)    are  then   constructed  for  the 
kl        k2 
remaining  pairs,    one  of  these    (again  the  one  associated  with  the  largest 

error)   necessarily  yields  the  best  a    . 

q 

Notice  that  the  construction  procedure  is  particularly  well- 
suited  to  parallel  computation.   The  best  unconstrained  approximations 
may  be  constructed  in  parallel,  following  which  the  prospective  a  's 
obtained  from  these  may  be  checked  for  all  of  the  other  curves  in 
parallel.   Finally,  if  several  pairwise  approximations  need  to  be 
computed,  they  may  also  be  done  in  parallel. 


12 


A  FORTRAN  program  implementing  the  algorithm  has  been  written 
with  the  aid  of  Dr.  Joseph  Garber.   This  program  consists  of  four  pieces: 

A.  A  main  program  which  does  the  bookkeeping  and  makes  the 
comparisons  described  in  the  first  paragraph  of  this  section. 

B.  A  subroutine  which,  given  an  exponential  factor  3,  finds  the 
best  approximation  aexp(3x)  to  a  single  function  g.  This  is 
the  subroutine  used  to  check  whether  an  a  obtained  from 

q. 

approximating  a  function   individually   does   in   fact   provide 
the  best   overall   exponential   factor. 

C.  A   subroutine  which   finds  the  best    exponential   approximation 

Sx 
ce       to   a  single   function   g. 

D.  A  subroutine  which   computes   the  best    simultaneous    approximation 
to   a  pair  of  functions,    given  that   the  characterization  is  by 
condition    (il)   of  Theorem  6. 

The   functions    {g   }    are  handled  in  tabular   form  throughout. 
This   is   of  some  importance   in   analyzing  experimental   data;    some 
approximation  programs    (e.g.    the  University  of   Illinois   subroutine 
library's  Remez   algorithm  program)    require   a  computable  g.      This 
necessitates  the  introduction  of  an   interpolatory  process    for  tabular 
data  and  tends  to   destroy  some  of  the  "minimax"   aspect  of  the 
approximation.      Also,   we   assume  that  the   functions    {g   }    are   essentially 
positive,    in  the   sense  that    for  each  k 

max  gk    (x)    -  min  g^   (x)   >  0. 
This    eliminates   the   complicating  possibility  that   the   zero   function  may 
appear  as   a  best   approximation  to  one  of  the   functions.      In  the   actual 
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coded  algorithm  and  some  of  the  details  to  follow  we  have  made  the  even 

more  restrictive  assumption  that  g,  >0  for  all  k.   However,  no 

theoretical  difficulties  should  arise  in  eliminating  this  restriction  so 

that,  for  example,  "noisy"  data  which  is  negative  at  a  few  points  may 

be  handled. 

The   simplest   of  Subroutines   B,    C,   D   is   B.      When   3   is   fixed, 

we  are   dealing  with  the  problem  of  finding  a  best   approximation  from 

the  one-parameter,    linear,    uni solvent   family    {aexp(gx):    as  R}.      Such 

approximations   exist,    are  unique,    and  are   characterized  by   a  two-point 

alternant.      That   is,    if  g  is  the   function  to  be   approximated,   then 

aexp(3x)   is  the  best   approximation   (on  the  given  interval)    if  and  only 

if  there  exist  two  points   x.  ,   x      in  the  interval   such  that    for  some 

constant   A 

g(x1)   -  aexp(3x]L)   =   A, 

(11)  g(x2)   =  aexp(6x2)    =  -A, 

and  p    E  g  -  aexp(gx)       |    =    |a|. 

The   single-exchange  Remez   algorithm,   which   is  known  to   converge   rapidly, 

is   used  for  Subroutine  B.      The  algorithm  proceeds   as    follows.      An 

initial   alternant   x      ,    x       is   guessed.      The  linear  equations    (ll)    are 

then  solved  for  a    and  A.      The  error    p  is  then  computed  and  compared 

with    |a|.      If  it   is   sufficiently   close,   the  procedure   is  terminated. 

Otherwise,    some  point  x*  at  which  the  maximum  error  is  taken  on  replaces 

either  x^      or  x       to   form  a  new  alternant   x.^    ,   x     ,    and  the  process   is 

repeated.      The  replacement   is  made  according  to  the  rule:      x*  replaces 

x     if  sgn   (g(x*)   -  aexp(3x*))   =  sgn   (g(x.)    -  aexp($x   ));   thus   the 
J  J  J 
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replacement   preserves   the   alternation  in   sign.      Convergence   is   very- 
rapid;    usually  only  one  or  two  iterations    are  required. 

Subroutine  C   is    somewhat   more   complicated  but    follows   along 

the  same  lines.      In  this    case  we  need  to   construct    a  best   approximation 

8x  2 

from  the  two-parameter,    nonlinear   family   {ae      :       (a, 3)    e  R   }. 

Characterization  is   now  by  three  points  of  alternation;   we   need  to    find 

points   x_<x  <x     such  that 

g(x    )    -  aexp(6x1)   =   X 
(12)  g(x    )    -  aexp(3x2)   =   -   A 

g(x    )    -  aexp(Bx   )   =    A 
and  PE||g  -  aexp(Bx)|      =    |a|. 

Subroutine  C   consists  of  a  single-exchange  Remez   algorithm,   but  there 
are  two  basic   difficulties.      First,    one  no  longer  has   a  convergence 
proof.      Second,   the  equations    (12)   to  be   solved  for  a  ,    6,    A   are  no 
longer  linear.      To     overcome   the  second  difficulty  we  follow  the 
procedures    suggested  by  Meinardus    [8,   p.    176   ff.].      By   selecting 
equally-spaced  points   for  the  initial   alternant    (i.e.   xp  -x_      = 
x_     -  xp     =   6 ) ,   we   find  that    addition  of  successive  pairs  of  equations 
from   (12)   to   eliminate   A   followed  by   division  of  the   resulting  two 
equations  to   eliminate  a  yields 

(13)  ■'-J441, 

gl        g2 

p 
where  we  have  used  g.  to  denote  g(x. ),  and  z  =  e  .   The  parameter  8 

is  immediately  obtained  by  taking  logarithms  in  (13),  and  values  of  a 

and  A  are  then  readily  computed.   Just  as  in  Subroutine  B,  the 

convergence  test  consists  of  checking  whether  |A|  is  sufficiently 
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close  to  p.      If  the  test   is   not   satisfied,    a  new  alternant   is  then 
formed  by  replacing  one  of  r     ,    x      ,    x       by  x*,    a  point  of  maximum 
deviation.      As   usual  the   alternation  of  signs  must  be  preserved;    for 
details   see  the  table  in   [8,   p.    107].      Equations    (12)   must   again  be 
solved  for  new  values   of  a,    3,    and   A,   but   now  the  points   of  the   alternant 
are  not    in  general   equally   spaced.      If  we  manipulate   Eqs.    (12)    just   as 
we  did  to  obtain   (13),   the  result   in  this    case   is   a  more   complicated, 
nonlinear  equation   in  3.      This   equation   can  be   solved  by  Newton's 
method,    following  which  a   and   A  may  then  be  obtained.      It   is   not,   however, 
really  necessary  to  solve  Eqs.    (12)   precisely  on   each  iteration  step. 
The  procedure   actually   coded  goes   as   follows.      Elimination  of  A   from 
successive  pairs  of  Eqs.    (12)   yields 

gl  +  g2  ~  a^exP(3x1)   +  exp(3x   )]  =  0 


(1>0 


l2  +  g3  "  atexP^x2)   +  exp(0x    )]   =  0. 


Supposing  that   the   parameters  a,    3   satisfying  these   are  reasonably   close 
to  the  a    ,    3      found  initially,   we  may   find  corrections    Aa  ,    A3   such  that 
(to   a  good  approximation)  a   =  a      +  Aa ,    3  =   3     +  A3  by  solving  the 
simultaneous   linear  equations 

(z1  +   Zg)    Aa    +  a°    (x^Z,    +   X2Z2>    A6   =   gl   +   g2   "  a°    ^Zl   +   Z2^ 
(z2  +   z    )    Aa    +  a°    (x^  +  X3ZJ   A^   =   gg  +   g     -  a°    (z      +   z    ), 
where   z.    =   exp(3  X.).      This    is,    of  course,    just   one   iteration   step   in 
the  Newton's  method  solution  of  the   system   (lU);   but   usually  this  one 
step   is    adequate   for  the   determination  of  a  new   alternant.       (incidentally, 
the   determinant   of   (15)    is   just 

a      ^z1z2^x2  "  xi^    +   Z2Z1^X3  "  X2^   +  Z1Z3    ^X3   "  xi^» 


16 


which   is   always   positive;   hence    (15)   may   always  be   solved.)      It    should 
"be  noted  that    since   Eqs.    (12)    are  not   satisfied  exactly,   the   convergence 
"test   in  this   case   does  not    involve    |x|    but   instead  the  quantity 

min   (|g(x.)    -  aexp(3x.  )  |  }  . 
i  =  1,2,3 
The   validity  of  this  test   comes    from  Theorem  85   in   [8],   or  alternately 
from  the  single- function  case  of  Theorem  1.      One  must,   of  course,    check 
that  the   error  signs    sgn    (g(x.  )   -  aexp(gx.))    alternate  on  x.^  ,    x   ,    x   . 
If  they  do  not,    additional   Newton's  method  cycles    (i.e.   the   solving  of 
(15)    for  successive   corrections)   will  be  necessary.      In  numerous  tests 
such   additional   cycles  were  never  found  necessary. 

Subroutine   C   as   a  whole  converged  extremely  rapidly,   requiring 
no   more  than  three  iterations    (i.e.    alternant   adjustments)  when  the 
given  functions    {g,  }  were  polynomials.      Convergence  was   just   slightly 
slower,    requiring  four  or   five  iterations,   when  the   (g,  }  were   exponentials 
with  random  noise   as   described  at   the   end  of  this    section. 

Two   approaches  were   implemented  for  the   construction  of  the 
best   simultaneous   approximation  to   a  pair  of  functions   g   ,    g    .      For 
reference  we   shall   call  the   resulting  subroutines   Dl   and  D2.      Subroutine 
Dl   is  based  on  the   fact  that    if   (3)   has    a  solution,    a  better  approximation 
may  be   constructed  from  that   solution.      That   is,  with   any  extremal    (x,   k) 
of  an  approximation   f(a,x)   there   is   associated  a  linear   inequality   (from 
(3))   of  the   form 

(16)      ak(x)    (3k  +  6     akx)   >  0. 
If  the  set   of  all   such  linear  inequalities    (associated  with  all   extremals) 
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has   a  solution  $,,$p>3    .  then   for  some   e>0   f(a',x)   provides   abetter 
approximation,   where   a'    =    (a      +   eg    ,  a      +   eB~,  a      +   eg    ) .      Thus  by 
iteratively   searching  for  extremals,   solving  inequalities    (l6),    and 
correcting  the  approximation,   one  may  hope  to  arrive   eventually   at 
a  best   approximation.      In  trials   this  method  has   never  failed  to 
converge,    although  even   for  simple  linear  approximation   convergence 
theorems   for  this   type  of  algorithm  do   not   seem  to  exist.      Convergence 
is    slow,   however,    and  the   rate  of  convergence   is   highly  dependent   upon 
the  method  of  choosing  the   adjustment   parameter   e.      Thus   a  complicated 
e  may   cut    down  on  the  number  of  iterations   required,   but  the  total  work 
involved  may  not   change  much.      Subroutine  Dl  was    abandoned  after  trials 
showed  that    it   required  something  like   four  times   as   many   iterations   and 
ten  times   as  much  time  as   Subroutine  D2.      It   is,  however,   worthwhile  to 
have  gained  some   experience  with  the  straightforward   "inequality"  method, 
since   it    is   easily   coded  and  may  work  in   some   situations    (e.g.    for  more 
complicated  approximating  families)    in  which   a  more  sophisticated 
algorithm  is   impossible. 

Subroutine  D2   is   a  Remez-type   algorithm,    very   similar  to 
Subroutines   B  and  C.      We  need  to   find  an   f  =    (cu    exp(Bx),  a    exp((3x)) 
such  that   at   four  points,    x^      <   x       and  x       <   x      ,  the   equations 

g1(x11)    -  a1  exp(6xi;L)   =    A 
g-L(x12)   -  a      exp(8x     )   =  -   A 


(IT) 


g2(x21)   -  a2   exp(6x21)   =   -  A 
g2(x22)    -a2   exp(gx22)   =  A 


hold  for  A  such  that  |a|  =  N(g  -  f).   Condition  (il)  of  Theorem  6  will 
then  be  satisfied.   Elimination  of  A  from  pairs  of  these  equations  results 
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in  the  three  equations 

(a)      g  +  g       =  a    [expCBx^)   +   exp(Bx^    )] 

(18)  .(b)      g21  +  g22  =  a2[exp(6x21)   +  exptgx^)] 

(c)      g±1  +  g21  =  a   exp(Bx^    )   +  a2exp(3x21), 

where  here  s. .  denotes  g.(x. t). 
Bij         Bi   ij 

If  (l8a)  and  (l8b)  are  solved  for  a  ,  a  ,  respectively,  and  these  values 
are  substituted  into  (l8c),  a  single  equation  in  3  is  obtained: 

£>11        "12  ^21        ^22 

(19)  H(3)    E  1   +  exp[6(x12  -  xi;l)]    +  1  +  exp[6(x22  -  x^)]    "  gll   "  g21  =  ° 

The  algorithm  then  proceeds   as    follows.      An  initial  set  of  points 
S  =    {x      ,    x      ,   x      ,    x22^   is    chosen   so  "that   x  _   -  x^      =  x_2  -  x?    .      Equation 
(19)    is   in  this   case   easily   solved  explicitly   for  the   initial   estimate  of 
3,    and  initial  values  of  a    ,   a    ,    A  are   subsequently  obtained  from   (l8)    and 
(IT).      The   usual   convergence  test — whether    |a|    is   sufficiently   close  to 
the  maximum  deviation  N(g  -   f)   of  the   approximation — is  then  made,    and  if 
it    fails   the  set   S  must  be   adjusted.      Equation    (19)    is  then  solved  for  a 
new   3  by  Newton's   method,    values  of  a    ,  a    ,    A   are  obtained  as   above,    and 
we   are  again   ready  to  test   for  convergence.      It   should  be  noted  that    (19) 
has   a  unique   solution,    since  H'(3)    <  0   and  H(3)  ■+■  ±    (g,,    +  g?1  )    as    3  ■+  +  °°. 

Adjustment   of  the  proposed  alternant  S  is  the  most  troublesome 
aspect   of  the  algorithm.      By  analogy  with  the  usual  procedure,   we  would 
like  the  new  S    (i)   to   contain   a  point  x*  of  maximum  deviation,   and   (ii) 
to  be  such  that  the  error  signs   alternate   in  the   desired  way   at   these 
points.      A  maximum  deviation  point   x*   is  therefore  located,    and  an 
attempt   is  made  to   replace  one  member  of  S  by  x*  in  such  a  way  that    (ii) 
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is   satisfied.      This   attempt  may  fail,    for  example   in  the   situation 
sketched  in  Fie.    1. 


We  see  that  when  we  replace  l      by  x*  to   get   the  proper  sign   alternation 
on  x*  <   x.     ,   the  required  sign   alternations  "between   curves    are  not 
satisfied.      To   satisfy  this   requirement,   multiple  replacements  must  be 
made.      Thus,    looking  at   Fig.    1,   we   find  that   either   {x*,   x      ,   y,*,    x     } 
or   (x*,    x      ,    x      ,  y  *}  would  be   satisfactory   as   the  new  set   S.      In   a 
slightly  more   complicated  case,   both  x        and  x       may  need  to  be  replaced, 
For  example,    if  g     -   f     is   as   shown   in  Fig.    1,  but   g     -   f     is   as    shown 
in  Fig.    2,    an   acceptable  new  set   S  would  be    {x*,   x., -,  ,   y-,*»  Y9*^' 
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Finally,  there  is  the  possibility  that  g  -  f  is  something  like  that 
shown  in  Fig.  3. 


Clearly  there  is  no  way  in  which  x*  may  be  inserted  into  the  set  S  so 
that  the  proper  sign  alternations  occur.   Looking  at  Fig.  1,  we  see 
that  a  small  shift  of  t.  to  the  left  may  have  the  proper  effect  of 
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evening  out  the  errors.      The  program  accordingly   carries  out   such   a 
small   shift,    after  first  printing  out   a  warning  message  that   an  alternant 
containing  x*   could  not  be  constructed.      This  possibility   causes   a 
serious   difficulty   in  attempts  to   analyze  the   algorithm  theoretically; 
it   is,   however,    not   very  likely  to   occur,    particularly   if  the   curves 
fitted  are   simply  "noisy"   exponentials.      In  33  trails    for  polynomial  g 
this  possibility  never  did  occur,    and  convergence  of  the   algorithm  was 
very   rapid,    requiring  on  the  average  U.l   iterations. 

Most   of  the  tests  of  the  program  were  carried  out  on   sets  of 
functions    (gv)  which   consisted  of  low-degree  polynomials  tabulated  at 
20   equally   spaced  points  on   [0,    l].      It  was    felt  that   such   functions 
would  adequately  test  the  main   features   of  the   algorithm  and  give 
reasonable   data  on  the  number  of  iterations   required  by  the  various 
subroutines,    etc.      The  program  as   a  whole   seems  to  be  quite  efficient; 
the  best   simultaneous   approximation  to  one   set  of  seven  polynomial 
functions  was  obtained  in  an  "execution  time"   of  under  1.5   seconds  on 
an  IBM  360/91.      During  the  course  of  this   computation,   Subroutine  Dl 
was    called  l6  times,    Subroutine  C    7  times,    and  Subroutine  B  ^7  times. 
Selected  excerpts   from  the  output   of  this  test    are  given  in  the  Appendix. 

Another  type  of  test   is  perhaps   of  more   interest    from  the 
point  of  view  of  practical   applications.      We  artificially  generated 
"experimental"   data  by   adding  random  errors   e    (|e|    <_. 01)   to   a  set   of 
three  exponential   curves  of  the     form  a e       with  a   =  0.5,   1.0,   1.5. 

("Data"   points  were   computed  for  20   equally   spaced  x-values  on   [0,2].) 

-h 
Our  program,   which   identified  extremals  with   a  tolerance  of  10      , 
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recovered  the   exponential   factor  £    =  -1   as   -1.0000.      The  traditional 
way  of  analyzing  exponential   data    (least -squares   straight -line  fit  to 
the  logarithms  of  the  function  values)   led  to  3   =  -0.9988.      The 
difference  is  largely   ascribable  to  the  weighting  induced  "by  taking 
logarithms;    direct   least-squares    fitting  of  exponentials   is,   however, 
a  troublesome  nonlinear  problem  even  for   a  single   function.      Simultaneous 
uniform  approximation   appears  to  be   a  very  practicable  alternative. 
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APPENDIX 
Excerpts    from  output   of  test:      simultaneous   exponential  fit  to  the 
functions 

gx  =  5  -   3x 

g2  =  h  -  3x 

g3  -  U  -  3x2 

g^  =  5  -  3x 

g5  =  7  -  2x 

gg  =  6  -   3x 

g~  =  7 
tabulated  at   20   equally-spaced  points  on    [0,   1]. 

Step  1.      Subroutine  C   is   used  to   find  the  best   exponential 

approximation  to  each  &     in  turn.      Below  are  two  examples   of  the  output 

of  Subroutine  C    (identified  as    "BEANR"   in  the  program).      Note  that   the 

point  of  maximum  deviation   (x*  or  XMAX)    and  the  points   of  the  proposed 

alternant    (Yl,    Y2,   Y3)    are  given  both   as  their  actual  values    and  as 

index  values    (IMAX,    II,    12,    13)   locating  them   in  the  tabular  data. 

The   numbers   SIGMAX  and  NEWSIG  actually   are  only  keeping  track  of  the 

error  signs    (at  XMAX   and  Yl,    respectively).      Parameters  a,    (B   are 

identified  as  A,    B  respectively.      Computations    are   for  g     and  g   . 

BERNR  SUBROUT I NE  CflLi  .ED 

INIT.IRL   VRLUE3   -  H:    9.  51411D  91        B:— S.  84y~>5D   99        LfiMBDR   ''--5IG"1    :  -9.  141060   99 

li:      i        Vi :    6.  88999D  66  12:    19        ¥2:    6.  473680   66  12:    19        Y3 :    6.  94737D   66 

ITERATION  NUMBER      1 

RHO  :  6.  198430  86     IMFltf:  26   KMRX :  6.  i9988D  81 

SIGMAX:  6.  i9843D  68   NEW  SIG:  6.  141860  66 

li:   i   Vi:  8.  868680  88     12:  18   ¥2 :  8.473680  86     13:  28   V3 :  8.188890  81 

A  :  8.  5iS14D  91   B  :  -8.  878320  88 


hCVlJ-Qc. 


'1  >  :  8  ±61  371 

F  v  V£ )  ■■  G  ( V2  '?  ""9  i  61  73[ 
F<V3)-Q<V3)  :  8.  161661 
L.  AM  BO  A      6   16 1330   89 


>  86 
88 


25 


I ! ERA ! I  ON  NUMBER  2 

RHO:  @.  161660  86     I  MAX:  26   KMflX:  9.  18800D  91 

SIGMFlM:  8.  161660  00   NEW  SIQ:  0.  14186D  ©0 

II:   1   VI:  9.  68888D  00     12:  19   V2 :  9.  47368D  99 

fl  :  8.  516i5D  91   B :  ~9.  97044D  80 

FCt'D-G^VI.')  :  9.  161460  89 

h  <  V2  >  -Q  (  V2  >  :  -9.  16146D  98 

FdV3?~CKV33  :  9.  1S14SD  99 

LAMBDA:  0.  16146D  88 


ITERATION  NUMBER  3 
RHO  :  8.  16146D  99 


;:mh:'-  2^   xMfiK  9  i08$gD  si 


13:  29 


V3: 


1998£ 


91 


END  OF  BEANR  SUBROUTINE 

BEFlNR  SUBROUTINE  CALLED 

INI  I IAL  VALUES  -  A;  8.  41958D  81 
II:   1   VI:  0.  88808D  88     12 


B: -8.  119410  81   LAMBDA  <-5IG; 

18 


r'2 :  o.  4r"j:6SD  88 


'K'MA!K-   9  16008D  81 
8.  195760  88 
12:  18   V2:  9.  47368D  00 


ITERA1 ION  NUMBER  i 

RHO:  8.  271200  88     IMAM:  28 
SIGriliX:  8.  271280  08   NEW  SIG 
II:   1   VI:  0.000O0O  OO 
A:  O.  422640  81   B:-8.  123670  k 
FCt'D-GCt'I?  :  O.  226390  OO 
h <V2>-G<V2>  :-8.  22626D  80 
F'::V3>-G<V3>  :  0.  227140  OO 
LAMBDA:  8.  226260  OO 


ITERATION  NUMBER  2 

RHO:  0.227140  OO     IMAX:  20   XMAX:  0.  108880  81 

SIGMAH:  6.  227140  88   NEW  SIG:  0.  19576D  "80 

II:   1   VI:  0.  OOOOOO  98     12:  10   V2 :  8.  4736S0  88 

A  :  8.  422660  81   B  :  -8.  1237i0  61 

F < VI > -G < VI >  :  6.  226640  88 

h  '.:V2>-G<V2>  :  -8.  226640  86 

F(V3>-G<V3>  :  6.  226640  86 

LAMBDA:  6.  226640  86 


ITERATION  NUMBER  3 

RHO  :  0.  22664D  88 

END  OF  BEANR  SUBROUTINE 


IMAM:  28 


MA .•■'  9.  18800D  Ol 


i: 


19 


20 


13 


-0.  195760  00 
V3 :  0,  947370  68 


8.  10008D  Ol 


O  10O0UO  Ol 


The   seven  best   individual   approximations    are  then  tabulated: 


I 

Ad) 

8(1) 

1 

8 

51615D 

61 

-8. 

S78440  88 

111 

8 

422660 

01 

-8. 

1237iD  81 

3 

8 

4*7090 

61 

-6. 

166810  61 

4 

8 

5i-'8950 

01 

-6. 

74*370  66 

5 

8 

764170 

61 

-8. 

334110  66 

6 

6 

612540 

81 

-6. 

b  .-2890  68 

7 

6 

786600 

61 

6. 

666880  86 

RHO(J) 

6.  161460  68 

6.  226640  66 

8.  578880  88 

6.  769480  66 
8.  416830-81 

8.  12538D  66 

0.  886880  88 


Step  2.      For  each   of  the  seven  exponential    factors   B(l) 
determined  above,    Subroutine   B    (identified  in  the  program  as    "BEAGR") 
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is   used  to   determine  the  best   approximation    (to   each  of  the  other  gk 

in  turn)   of  the   form  Aexp(B(l)x).      Output    for  B(l)    is   given  below.      The 

maximum  error  here  is   identified  as    "E"   in  the   BEAGR  output   and  "RHO" 

in  the   final  tabulated  output. 

1=     l 

BEflQR  SUBROUTINE  CALLED  WITH  B=-8.  87844D  60 

ITERATION  NUMBER  i 

t:  0.  475820  00   LAMBDA:-©.  47582D  00   A:  0.  35242D  01   VI:  0.  00000D  00 


END  OF  BEAGR  SUBROUTINE  -  A:  8.352420  01   E:  0.  47582D  00 

BEAGR  SUBROUTINE  CALLED  WITH  B=-0.  87044D  00 


V2:  0.  10000D  01 


iitKHfiUH  M'Ji'lb'tK   X 

E  :  0.  70311D  08   LAMBDA  :  -0.  67043D  SS   A  :  0.  355980  01   VI :  0.  31579D  00 

V2  ■  R  iFiPinnrc  ri 
ITERATION  NUMBER  2 

t :  8.  68277D  00   LAMBDA  :  -0.  6S277D  00   A :  0.  40184D  01   VI :  0.  42105D  00 

V2  :  0.  10000D  01 
END  OF  BEAGR  SUBROUTINE  -  A:  0.  40i84D  81   E:  8.  68277D  00 

BEAGR  SUBROUTINE  CALLED  WITH  B=-0.  87844D  00 


I  itRAl ION  NUMBER   1 

E  ;  0.  11812D  51   LAMBDA  :  -0.  5S545D  00   A  :  0.  8x3120  0i   VI 

ITERATION  NUMBER  2  >c"' 

t :  8.  898540  80   LAMBDA  :  -8.  82459D  88   A :  8.  58246D  81   VI 

V2 

VI 

V2 


0.  10000D   01 

8.  473680  80 

0.  08800D  88 

8.  578950  88 
8.  88888D  80 


ITERA1 ION  NUMBER  3 

E:  8.  870930  88   LAMBDA : -8.  878930  88   A:  8.  587890  81   VI 

END  OF  BEAGR  SUBROUTINE  -  A:  8.  5B785D  81    E:  0.  S7053D  00 

BEAGR  SUBROUTINE  CALLED  WITH  B=-8.  87844D  88 

ITERATION  NUMBER  1 

t :  8.  14588D  81   LAMBDA :  -8.  14588D  81   A  :  8.  84588D  81   VI :  8.  18688D  01 

V2:  8.  00088D  88 
END  OF  BtAGR  SUBROUTINE  -  A:  8.  845B8D  81   E:  8,  14588D  81 

BEAGR  SUBROUTINE  CALLED  WITH  B=-8.  87844D  80 


1  ItKHI  1UM    MUl'IfcStK       X 

E  :  8.  89847D  88   LAMBDA :  8.  11114D  88   A :  8.  68985D  81   VI :  0.  36842D  88 

V2 :  8.  18880D  01 
ITERATION  NUMBER  2 

E :  S.  43S48D  S5   LAMBDA :  8.  34333D  88   A :  8.  634350  81   VI :  8.  00808D  88 

V2 ;  8.  10880D  81 

ITERATION  NUMBER  3 

E :  8.  4i287D  88   LAMBDA :  8.  41247D  88   A :  8.  64125D  81   VI :  8.  088880  88 

V2:  0.  68421D  00 
ITERATION  NUMBER  4 

E  :  0.  41273D  88   LAMBDA :  8.  41273D  88   A :  8.  641270  01   VI :  0.  800000  88 

V2:  8.  73684D  88 
END  OF  BEAGR  SUBROUTINE  -  A:  y.  64127D  81   E:  5.  412730  00 

BEAGR  SUBROUTINE  CALLED  WITH  B=-0.  87044D  88 
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ITERATION  NUMBER  1 

E :  0,  28677D  01   LAMBDA :  -8.  28677D  8i   fl :  8.  98677D  81   VI :  Q.  10000D  91 

V2'  0  00000D  00 
END  OF  BEflQR  SUBROUTINE  -  R;  0  98677D  01   E:  0.  28677D  61 


J 

R<L  J> 

RHOa..  J) 

1 

8. 

51S15D 

01 

8.  16146D  08 

2 

0. 

33242P 

5i 

0 

475B2D  88 

T^ 

0, 

48184D 

01 

8 

68277D  88 

4 

8. 

58789D 

01 

8 

87893D  88 

5  . 

8. 

84588D 

01 

0 

14580D  81 

6 

8. 

64127D 

81 

8 

41273D  88 

7 

8. 

986770 

81 

0 

28677D  81 

Step.  3.   A  Boolean  matrix  M  =  (m  )  is  constructed  with 

J-  J 

m       =  T  if  RHO(l,   J)    <_  RHO   (i)   +   e   (a  small  tolerance)    and  m       =  F 
i  J  •  ij 

otherwise. 


THE  COMPARISON  MATRIX 

TFFFFFF 
(-  l  F  F  F  F  F 
TTTFFFF 
T  T  F  II-  r  F 
F  F  F  F  T  F  F 
F  F  F  h  F  T  F 
F   F   F   F   F   F   T 


Step,    h .      If  for  some  i  the   1       row  of  this   matrix  consists 
of  all   T's,    B(i)    is  the  best   overall   exponential   factor  and  the 
calcination  terminates.      This   did  not  occur  for  the   example   given  here, 

Step  5-      Using  Subroutine  D2    (identified  in  the  output   as 
"PAIRR",   the  program  now   computes  best    simultaneous   approximations   to 
pairs   {g   ,    g   }   such  that    (referring  to  matrix  M)   both  M       and  M        are 

p      q  pq  qp 

F.      Printout   of  results    for  three  such  pairs    is   given  below.      This    is 
largely  self-explanatory,    except   that    (J  =  1,    K  =  2)    indicates  that   it 
is    for  the   first   member    (g   )    of  the  pair  that   the   maximum  deviation   is 
taken  on,    and   (J  =   2,   K  =  l)    indicates  that    it   is    for  the   second 
member   (g   ). 
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TEP  5 


fa  » 9*1 


PAIRR  SUBROUTINE  CALLED  r  .  . 

INITIAL  VALUES  -  ALPHA1:  8.  53299D  61   ALPHA2:  8.  38487D  Si   BETH : -U.  1B.S5.SD  Wl 
RHO  :  8.  36388D  66    LAMBDA  .  6.  2651 3D  88 

ITERATION  NUMBER   1 

IMAX'  ':'0   "-'MAX •  8  18888D  61   J:  2   K:  1   bIG:  8.  .3638BD  08 
111-"  I"   Xll:  8.  52S32D-81     H2:  18   X12:  6.  S3474D  86 
121:   2   X21:  8.  5263:2D-6i     122 :  28   X22:  8.  i8888D  81_ 
pji_pj_ipij_  •  g  54i53D  81   RLFHR2  :  5.  37887D  81   t;t  i  A  :  -8.  188920  8i 
LAMBDA '  6.  27263D  66 
RH8:  6.  41535D  88 


ITERATluN  NUMBEK  2  _ 

IMAX:   1   XMAX:  6.  66666D  66  J:  1   K:  2   SIQ:  8.  415y5D  WW 

111-   1   Mil:  8.  88868D  66  112:  18   Xi2:  6.  89474D  86 

121-   2   X21:  6.  52632D-81  122:  26   X22:  8.  18888D  61 

ALPHAiT  6.  52S81D  61   ALPHA2  :  8.  37682D  61    BE  I  A  :  -8.  18713D  81 
LAMBDA:  6.  23868D  66 
RHO :  6.  41314D  88 

ITERATION  NUMBER  3 

IMAK'  12   XMAX:  6.  57335D  66   J:  1   K:  2   iIG:-8.  4iyl4D  88 
111:'  1   Xll:  6.  88888D  66     112:  12   X12:  8.  57895D  68 
77-1  •   2   >;:Z1 '  6.  52S32P- 8i     122:  -26   X22:  6.  iSSSSD  81 

ALPHA!  :  6.  53267D  61   ALPHA2  :  6.  371180  61   BETA  :  -6.  18286D  61 
LAMBDA '  6.  32668D  66 
RHO  :  8.  37333D  88 

ITERATION  NUMBER   4 

IMAX:   5   XMAX:  6.  21653D  66  J:  2   K:  1   SIQ: -8.  3799»i>  88 

111-   1   Xll:  6.  88888D  68  112:  12  X12:  6.  57895D  88 

121-   5   X2i:  6.  21653D  66  122:  26  X22:  8.  16666D  81 

ALPHA!  :  6.  53365D  61   ALPHA2  -  6.  37721D  61   BETA  :  -8.  18376D  61 
LAMBDA :_ 8.  33643D  66 
RHO:  6.  33374D  68 

ITERATION  NUMBER  5 

IMAX:   6   XMAX:  8.  26316D  86   J:  2   K:  1   SIQ :-8.  33974D  88 
111:   !   Xll:  8.  88888D  68     112:  12   X12:  8.  57895D  66 
121:   6   X21 :  6.  26316D  86     122:  26   X22:  8.  18688D  81 
ALPHA!  :  6.  53371D  61   ALPHA2  :  8.  37761D  61   BE1  A  :  -8.  18382D  81 
LAMBDA:  6.  3371.3D  86 
KHO  :  6.  33713D  66 

END  OF  PAIRR  SUBROUTINE 


hs  >  St} 


FRIRR  SUBROUTINE  CALLED 

INITIAL  VALUES  -  RLPHA1 :  5.  35676D  81   RLFHH2:  6.  54878D  61   BETR : -6.  83558D  68 
RHO:  6.  16344D  61   LfiMBDR  : -8.  252210  5R 

ITERATION  NUMBER   1 

IMAX:  12   XMAX:  8.  57535D  86   J:  2   K:  I   SIG:-8.  163440  81 
111:   2   Xll:  8.  526320-61     112:  18   X12 :  8.  83474D  88 
I2i:   2   X2i :  6.  52632D-8!     122:  12   X22 :  6.  57835D  86 
ALPHA1 :  6.  36321D  81   ALPHA2  :  8.  566880  81   BETA  :  -8.  62542D  88 
LAMBDA: -6.  47726D  66 
RHO  :  6.  18283D  81 

ITERATION  NUMBER  2 

IMAX:  28   XMAX:  6.  16868D  61   J:  2   K:  1   SIQ:  8.  18283D  81 

IMPOSSIBLE  TO  USE  XMAX 

111:   2   Xll:  8.  52632D-81     112:  18   X!2:  8.  83474D  86 

121:   2   X21:  8.  526320-81     122:  16   X22 :  8.  78347D  66 

ALPHA1:  8.  37689D  81   ALPHA2:  8.  55790D  61   BE l A: -6.  ,-22980  88 

LAMBDA:  -8.  37H8D  66 

RHO :  6.  82519D  86 
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ITERATION  NUMBER  3 

I  MAX:  25   Xf'iAX :  5.  155550  51 
lii:   2   Mil:  8.  52632D-81 
121:   2   X21:  8.  526320-81 
ALPHA1 :  6.  36010D  81   RLPHR2 : 
LAMBDA : -8.  54725D  00 
RHO :  0.  954850  55 


•-'  :  x   c. :  d.        :=•  j.  >j :  @.  BeiSlSu  yy 
112:  28   X12:  8.  166660  81 
122:  16   X22:  8.789470  88 
8.  57998D  81   BETR:-8.  84475D  88 


ITERATION  NUMBER  4 

I  MAX  :   8   XMAX  :  0.  363420  5u 
111:   3   Xli :  0.  36342D  00 
121:   2   K21:  0.526320-01 
RLPHR1 :  8.  41263D  01   RLRHR2 : 
LAMBDA : -8.  649440  88 
RHO  :  0.  93140D  88 


J:  1   K:  2   SIG: -5,  954850  08 
112:  28   X12 :  8.  188850  81 
122;  16   X22:  8.  789470  88 
8.  592830  81   BETA  :  -8.  91693D  88 


ITERATION  NUMBER  5 

I  MAX:  li   XMAX:  8.578550  88 
111:   S   Mil:  0,  36842D  88 
121:   2   X21 :  8.  526320-51 
ALPHA!:  8.  369950  81   ALPHA2 : 
LAMBDA:-©.  712690  86 
RHO:  0.  965220  88 


•J  :  ei   (<•. :  x   i>Iu:-5.  5314tf0  85 
112:  28   X12:  6.  186880  61 
icici:  lei   Xd'2 :  8.  578950  55 
6.  596520  81   BETA: -6.  822660  86 


ITERATION  NUMBER  6 

I MAX:   1   XMAX:  6.  686660  66 
Hi:   8   Mil:  6.368420  88 
121:   1   X21:  6.866880  66 
ALPHA! :  5.  378540  51   RLFHA2 : 
LAMBDA  :-0.  746250  66 
RHO  :  8.  746250  86 

END  OF  PRIRR  SUBROUTINE 


FAIRR  SUBROUTINE  CALLED 

INITIAL  VALUES  -  ALPHA1:  8.318560  81 
RHO  :  8.  1^.7540  51   LAMBOA  :  -8.  855280  88 


J:  2   K:  i   SIG:  8.965220  88 
112:  28   X12:  8.  186860  61 
122:  12   X22:  6.  578950  66 
5.  574620  51   BETA  :  -5.  773680  w 


[93>9r] 


ALPHA2:  8.797670  81   BETA : -8.  291520  80 


ITERAl ION  NUMBER  i 

I MAX:  26   XMAX:  6.  i68880  81 
111:   2   Xli:  6.  526320-6i 
121:   2   X2i :  6.  526.s20~6i 

ALPHA! :  5.  258520  5! 
LAMBOA  :  -5.  156820  5! 
RHO:  6.  129260  81 


J:  1   K:  2   SIG:  6.  137940  81 

1 12  :  26   X12  :  8.  166660  61 


122:    18 


ML^MHii  :     y.   St'^yLJ    ox 


;94740   86 


bt  i  i-t :  — &  st-OtZfu  otj 


ITERATION  NUMBER   2 

I MAX:  26   XMAX :  8.166680  81 
111:   2   Xli:  6.526320-61 
121:   2   X21 :  5.  526320-5! 
ALPHA!:  6.  294840  01   ALPHA2 
LAMBOA  :-8.  11027D  81 
RHO  :  6.  124760  01 


J:  2   K:  1   SIG: -6.  129260  81 
112:  28   X12:  8.168880  81 
122:  26   X22 :  5.  155550  51 
8.  824760  01   BETA: -6.  335350  86 


ITERATION  NUMBER   3 

I MAX:   1   XMAX:  8.868660  86 
111:   2   Xli:  8.  526320-61 
121:   i   XS1:  0,  tiOOfOO  55 
ALPHA! :  8.  292250  51   ALFHA2 : 
LAriBOA  :  -5.  111790  81 
RHO  :  6.  114730  81 


J:  2   K: 
I 12 :  28 

1 22 :  25   > 
":1!750  51 


0. 


i   SIG:  8.  124780  81 
X12 :  6.  188880  61 
X22 :  6.  186550  61 

BETA  :  -6.  322140  56 


ITERAl ION  NUMBER  4 

I  MAX:   4   XMAX:  6.  is  .''890  86 
Hi:   4   Xli:  6.  i57890  66 
I2i;   1   X21 :  6.  668660  66 
ALPHAi :  6.  294480  61   ALPHA2 : 
LAMBOA: -6.  11277D  61 
RHO:  6.  112770  61 


J:  1   K:  2   bIG:-6.  114730  01 
112:  26   X12:  6.  188680  61 
122:  28   X22:  6.186660  61 
8.  81277D  61   BETA: -6.  325820  68 


END  OF  PAIRR  SUBROUTINE 
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After  all  pairwise  approximations  have  been  computed,    a 
summary  of  the  results   are  printed  out   and  the  pair  yielding  the 
largest   p    =  N(g  -   f)   is   identified.      In  this   example,   the   characterizing 
pair  is    {gy    g   }. 


I 

J 

RLPHRKL 

•J> 

RLPHR2<I. 

J) 

BETfKLJ) 

RHCKL  J 

> 

i 

2 

0. 

53371D 

8i 

0 

37761D 

01 

"0 

10382D 

01 

8. 

33713D 

08 

i 

£i 

9. 

442i±D 

01 

0 

75789D 

01 

-0 

53900D 

00 

0. 

D789tiD 

88 

l 

6 

0. 

48432D 

01 

0 

62515D 

01 

-8 

76597D 

00 

8. 

25i51D 

08 

1 

7 

0. 

40088D 

01 

0 

80900D 

01 

-0 

28768D 

00 

0. 

10000D 

81 

2 

3 

0. 

32353D 

01 

0 

77647D 

01 
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Finally,  Subroutine  B  is  used  to  find  the  best  coefficients 
C   in  the  approximation  of  the  g,'s  by  C   exp  (3x),  where  the  B  is  as 
obtained  above.   (in  this  example,  6  =  BETA  (3,  T)  =  -0.32502.) 
Errors  p,  of  these  approximations  are  also  printed  out. 
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